package com.snap.vseries.analyze;

import android.util.Log;
import android.util.SparseArray;

import com.orhanobut.logger.Logger;
import com.snap.vseries.bean.AnalysisProtocol;
import com.snap.vseries.bean.ChannelBean;

import org.apache.commons.math3.exception.MaxCountExceededException;
import org.apache.commons.math3.fitting.PolynomialCurveFitter;
import org.apache.commons.math3.fitting.WeightedObservedPoints;
import org.apache.commons.math3.fitting.leastsquares.ParameterValidator;
import org.apache.commons.math3.linear.RealVector;
import org.apache.commons.math3.ml.clustering.CentroidCluster;
import org.apache.commons.math3.ml.clustering.DoublePoint;
import org.apache.commons.math3.ml.clustering.KMeansPlusPlusClusterer;
import org.apache.commons.math3.ml.clustering.MultiKMeansPlusPlusClusterer;

import java.util.ArrayList;
import java.util.Arrays;
import java.util.HashSet;
import java.util.List;
import java.util.Set;

public class MoleAlgo {


    private static final String TAG = "MoleAlgo";
    private static final int LINE_FIT_PRE_EXCLUDE = 3;
    private static final float NORMAL_AFI_BASE_RATE = 0.12f;//afi和基数的比例关系，2000基数为例  afi临界值360；
    private static final float END_POINT_RATE = 1.05f;//afi和基数的比例关系，2000基数为例  afi临界值360；
    private static final float SLOPE_BASE_RATE = 0.002f;//斜率和基数的比例关系，2000基数为例 直线拟合最大斜率为4；

    public static void main(String[] args) {
        double[] param = {16362.113671875002, 133.42715938541082, 5838.193708357465, 1.0, 22.560825803447607, 0.21902525863301764, 49.18808549607251};
//        double[] param = {1.4707*20, 0.0093, 10.9421*20, 0.1081, 49.1868, 0.7858, 35.9089};
//        double[] param = {5000, 5.6, 2000/5, 0.1, 23, 0.78, 32};
        CurveFunction a = new CurveFunction(param);
        int maxIndex = 100;
        double max = 0.0;
        for (int k = 100; k < 4000; k++) {
            double current = a.value(k / 100.0);
            if (current > max) {
                max = current;
                maxIndex = k;
            }
        }
        System.out.println(Arrays.toString(param));
        System.out.println(maxIndex);
        System.out.println(max);
        System.out.println(1 / max);

        //扩增高倍率  以这个为标准  ct延后了   0.4 有点影响了
        //[1.4, 5.6, 20.0, 0.1, 23.0, 0.78, 32.0]
        //2858
        //0.27279432586930724
        //3.665765396011532


        //扩增高倍率 c=c/a;
        //[1.4, 5.6, 14.285714285714286, 0.1, 23.0, 0.78, 32.0]
        //2898
        //0.2571782842610944
        //3.8883531822024784


        //扩增低倍率   以这个为标准  c除了后 ct延后了不多 0.08
        //[1.4, 5.6, 3.0, 0.1, 23.0, 0.78, 32.0]
        //3027
        //0.12104532051185808
        //8.261368516943502


        //扩增低倍率 c=c/a;
        //[1.4, 5.6, 2.142857142857143, 0.1, 23.0, 0.78, 32.0]
        //3035
        //0.09119076459452252
        //10.966022759502843


        //[500.0, 5.6, 3000.0, 0.1, 23.0, 0.78, 32.0]
        //2245
        //0.3184988357156313
        //3.1397289027858193

//        [500.0, 5.6, 6.0, 0.1, 23.0, 0.78, 32.0]
//        2988
//        0.19119225065976364
//        5.230337508707668

    }

    public static void execute(SparseArray<ChannelBean> channelBeanSparseArray, List<AnalysisProtocol> protocols) {
        //按通道归类
        List<List<ChannelBean>> list = new ArrayList<>();
        list.add(new ArrayList<>());
        list.add(new ArrayList<>());
        list.add(new ArrayList<>());
        list.add(new ArrayList<>());
        for (int i = 0; i < channelBeanSparseArray.size(); i++) {
            ChannelBean bean = channelBeanSparseArray.valueAt(i);
            list.get(bean.getChannelIndex()).add(bean);
        }

        for (int i = 0; i < list.size(); i++) {
            List<ChannelBean> channelList = list.get(i);
            for (int j = 0; j < channelList.size(); j++) {
                single(channelList.get(j), protocols.get(i).isAutoThreshold(), protocols.get(i).getThreshold());
            }
        }
    }

    private static float single(ChannelBean bean, boolean autoThreshold, float threshold) {
        if (bean.getOrigin() != null && bean.getOrigin().length >= 10) {
            if (lineFit(bean)) {
                //直线拟合成功
            } else if (curveAmpFit(bean, autoThreshold, threshold)) {
                //曲线拟合成功
            }
        }
        return 0;
    }

    private static boolean lineFit(ChannelBean bean) {
        WeightedObservedPoints points = new WeightedObservedPoints();
        for (int x = LINE_FIT_PRE_EXCLUDE; x < bean.getOrigin().length; x++) {
            points.add(x + 1, bean.getOrigin()[x]);
        }
        PolynomialCurveFitter fitter = PolynomialCurveFitter.create(1);
        try {
            double[] params = fitter.fit(points.toList());
            if (params[1] < params[0] * SLOPE_BASE_RATE) {
                Log.i(TAG, "判定为直线啦");
                return true;
            } else if (params[1] * bean.getOrigin().length < 100) {
                Log.i(TAG, "判定为直线啦");
                return true;
            } else {
                return false;
            }
        } catch (Exception e) {
            return false;
        }
    }

    private static boolean curveAmpFit(ChannelBean bean, boolean autoThreshold, double threshold) {
        float[] copy = bean.getOrigin().clone();
        double[] initParams = getInitGuess(copy, bean);

        final MyCurveFitter fitter = new MyCurveFitter(
                new DoubleSigmoid.Parametric(),
                initParams, 30000,//拟合最大迭代次数,设置多少对总耗时影响不大，// 正常的数据在100次内都能迭代出来，异常数据 10000次都可能失败
                new MyParameterValidator());

        //1.前10个的权重降低
        //2.权重逐渐递增
        //3.前n个数据不参与拟合
        WeightedObservedPoints points = new WeightedObservedPoints();
        for (int i = LINE_FIT_PRE_EXCLUDE; i < copy.length; i++) {
            points.add(i + 1, copy[i]);
        }

        try {
            double[] curveFitParams = fitter.fit(points.toList());

            if (curveFitParams != null) {
                //给参数 de 和 fg 排顺序(对计算Ct 求导等不产生任何影响）
                if (curveFitParams[4] > curveFitParams[6]) {
                    double temp = curveFitParams[4];
                    curveFitParams[4] = curveFitParams[6];
                    curveFitParams[6] = temp;
                    temp = curveFitParams[3];
                    curveFitParams[3] = curveFitParams[5];
                    curveFitParams[5] = temp;
                }
                Log.e(TAG, "曲线拟合结果：" + Arrays.toString(curveFitParams));
            }

            if (curveFitParams == null) {
                Log.i(TAG, "曲线拟合失败");
                return false;
            } else if (curveFitParams[2] < curveFitParams[0] * NORMAL_AFI_BASE_RATE) {
                Log.i(TAG, "afi过低，拟合不合格：" + curveFitParams[2]);
                return false;
            } else if (curveFitParams[3] + curveFitParams[5] < 0.15) {
                Log.i(TAG, "两个锐度都过于小，可能是斜线，拟合不合格");
                return false;
            } else {
                //计算ct值
                CurveFunction curveFunction = new CurveFunction(curveFitParams);
                int maxIndex = 100;
                double max = 0.0;
                for (int k = 100; k < bean.getOrigin().length * 100; k++) {
                    double current = curveFunction.value(k / 100.0);
//                    current = Math.floor(current * 1000) / 1000;//会影响一定的精确性 每条线提前的距离不一样。影响cv、标准曲线等
                    if (current > max) {
                        max = current;
                        maxIndex = k;
                    }
                }
                float ctValue = maxIndex / 100f;
                Log.e(TAG, "ct值: " + ctValue);
                if (ctValue < 12) {
                    Log.i(TAG, "ct值过低，拟合不合格");
                    return false;
                }
                DoubleSigmoidMulti f = new DoubleSigmoidMulti(curveFitParams);
                if (f.value(bean.getOrigin().length) < END_POINT_RATE) {
                    Log.e(TAG, "最大循环处的倍率：" + f.value(bean.getOrigin().length));
                    Log.i(TAG, "倍数趋势不明显、、、、从尾部判断、、、、，拟合不合格");
                    return false;
                }
                Log.e(TAG, "曲线拟合成功");


                bean.setFitParams(curveFitParams);
                bean.setAmpFlag(true);
                //填充相对荧光数据
                float[] fit = new float[bean.getOrigin().length];
                for (int k = 0; k < fit.length; k++) {
                    fit[k] = (float) f.value(k + 1);
                }
                bean.setFit(fit);
                if (autoThreshold) {
                    bean.setCtValue(ctValue);
                } else {
                    bean.setCtValue(calcuManualCt(bean, threshold));
                }
                return true;
            }
        } catch (Exception e) {
            if (e instanceof MaxCountExceededException) {
                Logger.e("曲线拟合时超过最大迭代次数！");
            } else {
                e.printStackTrace();
            }
            return false;
        }
    }

    private static float calcuManualCt(ChannelBean bean, double thresholdValue) {
        if (bean == null || bean.getFitParams() == null || bean.getOrigin() == null) {
            return 0;
        }
        DoubleSigmoidMulti f = new DoubleSigmoidMulti(bean.getFitParams());
        for (int k = 100; k < bean.getOrigin().length * 100; k++) {
            double value = f.value(k / 100.0);
            value = Math.floor(value * 1000) / 1000;
            if (value >= thresholdValue) {
                return k / 100f;
            }
        }
        return 0;
    }

    /**
     * 初始猜测
     */
    private static double[] getInitGuess(float[] data, ChannelBean bean) {

        final double start;//起点很关键，因为回归算法对这个值不做任何修正
        int baseLineStart;
        int baseLineStop;
        if (data.length < 10) {
            baseLineStop = data.length;
        } else {
            Set<Float> set = new HashSet<>();
            for (float datum : data) set.add(datum);
            if (set.size() < 4) {
                baseLineStop = Math.min(data.length, 15);
            } else {
                baseLineStop = getBaseLineEnd(data);
            }
        }
        if (baseLineStop > 3) baseLineStart = 3;
        else baseLineStart = 0;
        bean.setBlStart(baseLineStart);
        bean.setBlStop(baseLineStop);
        double[] baseLineFitParams = baseLineFit(data, baseLineStart, baseLineStop);
        if (baseLineFitParams == null) {
            start = (data[4] + data[5] + data[6]) / 3;
        } else {
            start = baseLineFitParams[0];
        }
        final double end3Avg = (data[data.length - 1] + data[data.length - 2] + data[data.length - 3]) / 3;
        final double afi = end3Avg - start;

        double sum = 0.0;
        for (float datum : data) {
            sum += datum;
        }
        final double avg = sum / data.length;
        //最初版本
//        int firstAbove = 0;
//        while (data[firstAbove++] < avg) ;
        //枪毙原因：在前10个循环内有一个特别高的异常值时，firstAbove会取这个异常值，这种初始猜测会严重影响拟合效果

        //改进版本
//        int firstAbove = 10;
//        if (data.length > 10) {
//            while (data[firstAbove] < avg) {
//                firstAbove++;
//            }
//        }
        //发现了个数组越界问题，再修复后得到目前的
        int firstAbove = 10;
        while (firstAbove < data.length && data[firstAbove] < avg) {
            firstAbove++;
        }

        return new double[]{start, baseLineFitParams[1], afi, 0.1, firstAbove, 0.7, firstAbove};
    }

    /**
     * 聚类算法猜测基线截止点 且只是猜个大概。
     */
    private static int getBaseLineEnd(float[] data) {
        List<DoublePoint> dataList = new ArrayList<>();
        for (int i = 0; i < data.length; i++) {
            dataList.add(new DoublePoint(new double[]{data[i]}));
        }
        MultiKMeansPlusPlusClusterer<DoublePoint> kpp = new MultiKMeansPlusPlusClusterer<>(new KMeansPlusPlusClusterer<>(4), 10);
        List<CentroidCluster<DoublePoint>> clusterList = kpp.cluster(dataList);
        CentroidCluster<DoublePoint>[] clusterArr = clusterList.toArray(new CentroidCluster[0]);
        Arrays.sort(clusterArr, (o1, o2) -> Double.compare(o1.getCenter().getPoint()[0], o2.getCenter().getPoint()[0]));
        int[] sizeArr = new int[4];
        sizeArr[0] = clusterArr[0].getPoints().size();
        sizeArr[1] = clusterArr[1].getPoints().size();
        sizeArr[2] = clusterArr[2].getPoints().size();
        sizeArr[3] = clusterArr[3].getPoints().size();
        int result;
        if (sizeArr[0] > 5) {
            result = sizeArr[0];
        } else if (sizeArr[1] > 5) {
            result = sizeArr[0] + sizeArr[1];
        } else if (sizeArr[2] > 5) {
            result = sizeArr[0] + sizeArr[1] + sizeArr[2];
        } else {
            result = data.length;
        }
        return result - 5;
    }

    /**
     * 基线的直线拟合
     */
    private static double[] baseLineFit(float[] data, int baseLineStart, int baseLineEnd) {
        WeightedObservedPoints points = new WeightedObservedPoints();
        for (int x = baseLineStart; x <= baseLineEnd; x++) {
            points.add(x + 1, data[x]);
        }
        PolynomialCurveFitter fitter = PolynomialCurveFitter.create(1);
        try {
            return fitter.fit(points.toList());
        } catch (Exception e) {
            return null;
        }
    }

    public static class MyParameterValidator implements ParameterValidator {

        //参数修正器中的修正范围
        final int minLimit = 0;
        final int maxLimit = 50;

        @Override
        public RealVector validate(RealVector params) {
            double p0 = params.getEntry(0);
//            double p1 = params.getEntry(1);
            double p2 = params.getEntry(2);
            // afi: 0~5倍
            if (p2 > 5 * p0) {
                params.setEntry(2, 5 * p0);
            } else if (p2 < 0) {
                params.setEntry(2, 0);
            }

            //这个修正是错误的。会导致拟合的线和 原始的线完全不一致
//            //基线坡度:
//            if (p1 > p0 * SLOPE_BASE_RATE) {
//                params.setEntry(1, SLOPE_BASE_RATE * p0);
//            }

            double p3 = params.getEntry(3);
            if (p3 > 1.0) {
                params.setEntry(3, 1.0);
            } else if (p3 < 0) {
                params.setEntry(3, 0.1);
            }

            double p5 = params.getEntry(5);
            if (p5 > 1.0) {
                params.setEntry(5, 1.0);
            } else if (p5 < 0) {
                params.setEntry(5, 0.1);
            }

            double p4 = params.getEntry(4);
            if (p4 > maxLimit) {
                params.setEntry(4, maxLimit);
            } else if (p4 < minLimit) {
                params.setEntry(4, minLimit);
            }

            double p6 = params.getEntry(6);
            if (p6 > maxLimit) {
                params.setEntry(6, maxLimit);
            } else if (p6 < minLimit) {
                params.setEntry(6, minLimit);
            }
            return params;
        }
    }
}
